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Abstract 

We present a study of the self consistent Ornstein-Zernike approximation (SCOZA) for square- 
well (SW) potentials of narrow width 8. The main purpose of this investigation is to elucidate 
whether in the limit 5 — > 0, the SCOZA predicts a finite value for the second virial coefficient at 
the critical temperature B2(T C ), and whether this theory can lead to an improvement of the approx- 
imate Percus-Yevick solution of the sticky hard-sphere (SHS) model due to Baxter [R. J. Baxter, 
J. Chem. Phys. 49, 2770 (1968)]. For SW of non vanishing 5, the difficulties due to the in- 
fluence of the boundary condition at high density already encountered in an earlier investigation 
[E. Scholl-Paschinger, A. L. Benavides, and R. Castaneda-Priego, J. Chem. Phys. 123, 234513 
(2005)] prevented us from obtaining reliable results for 5 < 0.1. In the sticky limit this difficulty 
can be circumvented, but then the SCOZA fails to predict a liquid-vapor transition. The picture 
that emerges from this study is that for 6 — > 0, the SCOZA does not fulfill the expected prediction 
of a constant B 2 (T C ) [M. G. Noro and D. Frenkel, J. Chem. Phys. 113, 2941 (2000)], and that for 
thermodynamic consistency to be usefully exploited in this regime, one should probably go beyond 
the Ornstein-Zernike ansatz. 
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I. INTRODUCTION 



Narrow attractive interactions are ubiquitous in the realm of colloidal dispersions. A most 
relevant instance is represented by depletion-induced forces in mixtures of highly asymmetric 
hard spheres, a field that has received a fundamental contribution by Bob Evans 1|. Al- 
though these investigations have shown that the structure of effective depletion interactions 
is in general more complex than just a single attractive tail, a narrow attractive contribution 
at short distance is still to be expected in most cases, and its relevance to the modelization 
of both the equilibrium and non equilibrium phase behavior of colloidal systems has long 
been acknowledged. 

In view of this, it is somewhat disappointing that the difficulties experienced by most 
liquid-state theories when phase separation occurs, the foremost of which is lack of conver- 
gence in a (not necessarily small) neighborhood of the critical region, only become worse 
when the range of the attractive potential decreases to a small fraction of the particle size. 
However, it may be argued that this is not a fundamental difficulty, inasmuch as this class 
of interactions can be dealt with by the sticky hard-sphere (SHS) model [2|. This is ob- 
tained from the prototypical attractive interaction, i.e., the hard-core plus square- well (SW) 
potential t>sw( r ) given by 
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by taking the limits 5 — > 0, e — > 00 so that the second virial coefficient B 2 stays finite and 
non vanishing: 
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The stickiness parameter r is related to the second virial coefficient of the SHS model i?| HS 
by the relation 

73HS 

r= (3) 
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where I?? s = 2ira 3 /3 is the second virial coefficient of the hard-sphere fluid. As shown 

n 

by Baxter [2|, this model be can solved exactly within the Percus-Yevick (PY) integral 



equation Because of the way the limit ([2]) is taken, the solution does not depend on 
the temperature, but rather on the stickiness r or, equivalently, on I?f HS . The connection 
between the SHS model and an actual fluid at a certain temperature T can be made by a 
generalized principle of corresponding states that, although not rigorous, has been found to 
work remarkably well. According to this principle, usually referred to as Noro-Frenkel (NF) 
scaling J4], a pair potential v(r) consisting of a hard-core repulsion with diameter a and an 
attractive tail w(r) can be mapped into the SW potential (pQ), provided that the width 8 of 
the attractive well is determined so that its second virial coefficient -f^CO will be the same 
as that of the original interaction at the same temperature T, assuming that in both cases T 
is measured in units of the well depth. On the other hand, when 8 is small, by following the 
same procedure the SW potential can in turn be mapped into the SHS model. In this way, 
NF scaling can be combined with Baxter analytical solution of the SHS model to provide 
a unified description of fluids with narrow attractive interactions. For instance, the critical 
temperature T c of such a fluid can be evaluated by that at which its B2(T) is equal to the 
critical value of I?f HS . 

Although the SHS model has a perfectly regular solution within the PY approximation 
and has also been the subject of intensive simulation studies , the singularity embedded 
into its very definition should instigate a mild feeling of suspicion. In fact, twenty-three 
years after the appearance of the model, Stell j^] discovered that the latter from a rigorous 
statistical mechanical point of view is thermodynamically unstable whenever the number of 
particles is greater than eleven. However, the state of affairs is less serious than it might 
seem: the instability arises from a number of rather fancy particle configurations, that are 
readily destroyed by the slightest amount of polidispersity in the system. Therefore, the 
SHS model and its approximate analytical solution via the PY equation still remain a 
valuable tool for the interpretation of the phase diagram of colloidal dispersions subject to 
short-range interactions. 

There is, however, a significant limitation of this solution that negatively affects its pre- 
dictive power, i.e., its lack of thermodynamic consistency. This is particularly evident when 
the liquid-vapor transition of the model is considered: the coexistence curves determined by 
the compressibility and internal energy routes described in Sec. [Til differ markedly from each 
other, the discrepancy between the critical densities amounting to nearly a factor three, and 
neither of them agrees quantitatively with the simulation results In this respect, it is 
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also worth recalling that the PY approximation upon which the analytical solution of the 
SHS is based, is expected to become less and less accurate as the density is increased, while 
on the other hand the critical density increases as the range of the interaction decreases, 

nn 

and for 8 — > the limiting value predicted by the simulations p, [10[ is considerably larger 
than that expected for Lennard- Jones like fluids. 

Recently, a different approach has emerged, which is based upon enforcing consistency be- 
tween the aforementioned energy and compressibility routes to the thermodynamics, namely, 
the self-consistent Ornstein-Zernike approximation (SCOZA) llNl5|. This has been applied 
especially to the case in which v(r) consists of a hard core and a Yukawa (HCY) attractive 
tail: 

foo r < a , 

^HCY(r) = { a (4) 
-e-e z(r/a 1J r > a , 



where z is the inverse-range parameter of the interaction. For the HCY potential, the SCOZA 
has shown a remarkable robustness even for very short-range tails, while still remaining 
accurate [ijj] . On the other hand, one might argue that accuracy cannot be expected 
to persist for attractions of arbitrarily short range, since after all, as recalled in Sec. HIl 
the SCOZA assumes that the direct correlation function c(r) of the fluid is linear in the 
attractive part w(r) of the interaction, and its low-density expansion does not even give 
back the correct expression of B 2 (T). In 16] it was suggested that for narrow attractive 
interactions it could be more appropriate to have a c(r) which depends linearly not on w(r), 
but rather on its Mayer function e~^ w ^ — 1, where = l/{k-oT) is the inverse temperature. 
If this prescription were adopted in the SCOZA, one would obtain a "non-linear SCOZA" 
bound to become exact at low density. 

Now, the SW potential is quite special in this respect, because both w(r) and its Mayer 
function are constant in the well region and vanish elsewhere, and therefore they have the 
same functional form as far as they dependence on r is concerned, the only difference being 
their (state-dependent) amplitude. On the other hand, in the SCOZA this amplitude is 
determined by an exact thermodynamic consistency condition. As a consequence, for the 
SW potential the usual SCOZA must coincide with its non-linear counterpart, and hence 
predict the exact B 2 (T) at low density, unlike for a generic tail potential. It is then tempting 
to investigate how the SCOZA will behave for narrow SW potentials: will the quite special 
status of this potential enable it to recover a finite B 2 (T C ) in the limit 5 — > 0? This possibility 
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is quite enticing, because it would allow one to fix the lack of thermodynamic consistency 
that is the main defect of the PY solution of the SHS fluid (2L possibly obtaining a better 
agreement with the simulation results for the phase diagram [7|. 

With this in mind, we have performed a study of the SCOZA for narrow SW interactions. 
Such a task requires a fully numerical solution of the SCOZA, because one cannot take 
advantage of the wealth of analytical results available for the HCY potentials \y\ - In 
fact, this problem was already tackled in jisj], where the numerical solution of the SCOZA 
for a generic tail interaction w(r) was worked out and applied to the SW potential. The 
investigation performed there showed that on narrowing the well, the results become more 
and more dependent on the choice of the upper boundary p of the density interval where 
the SCOZA is solved, at which a high-density boundary condition needs to be specified. In 
order to get rid of such a dependence, one needs to move po to higher and higher values, 
but eventually this is prevented by a lack of convergence of the numerical algorithm. On 
the basis of the conclusions of that work, we have then reconsidered the numerical solution 
of the SCOZA paying special attention to the problem of obtaining convergence at high 
density, and developed a numerical algorithm that remains convergent up to po — l-4cr~ 3 , 
which is near the density at close packing. By this, we have managed to obtain results for 
wells down to 5 = 0.1, which however is not small enough to represent reliably the behavior 
of the SCOZA in the limit 5^0. 

We have then resorted to a different strategy, namely, we have taken the sticky limit ([2]) 
of the SW potential from the outset, and imposed thermodynamic consistency on the top 
of this. Again, the results turn out to be strongly dependent on the high-density boundary 
condition, but in this case the simpler form of the SCOZA equation allows one to circumvent 
this problem by resorting to a solution technique which requires only the behavior at low 
density to be specified. Unfortunately, according to this solution the SCOZA fails to have a 
liquid- vapor transition in the sticky limit (|2]). 

However disappointing, the full picture appears to be consistent enough to indicate that 
the SCOZA critical temperature will decrease too rapidly to give a finite value of B2(T C ) for 
5—7-0. On the other hand, the approximation which is introduced in the SCOZA, namely, 
that the contribution to the direct correlation function c(r) due to the tail interaction w(r) 
has the same range as w(r) itself, and therefore vanishes wherever w(r) does; is made 
also in the PY equation, which moreover is thermodynamically inconsistent. This suggests 
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that using thermodynamic consistency as a way to improve the PY solution for the SHS 
potential fl2]) would require to give up the ansatz on the range of c(r). 

The plan of the paper is as follows: in Sec. [TTJ the algorithm employed for the fully numer- 
ical solution in the case of a generic attractive tail w(r) is described; in Sec. IHII our results for 
SW potentials of several amplitudes 6 are presented, and compared with numerical simula- 
tions; the sticky limit of the SCOZA is considered in Sec.lIV( our conclusions are presented in 
Sec. |VJ At the end of the paper, three Appendixes deal with some technical details pertinent 
to Sec. [TTJ Appendix IA1 presents the finite- difference algorithm used to integrate the SCOZA 
partial differential equation (PDE); Appendix iBl concerns the minimization algorithm which 
we used to optimize c(r) inside the repulsive core of the interaction; Appendix O describes 
the technique adopted for the treatment of the hard-sphere part of the correlations at high 
density. 



II. THE SCOZA FOR A GENERIC ATTRACTIVE TAIL POTENTIAL 



We are interested in a simple model fluid of particles interacting via a two-body, spheri- 
cally symmetric potential v(r) that consists of a hard core with diameter a and an attractive 
tail w(r): 

+00 r < a , 

w[r) r > a . 

As is customary in integral-equation theories, the SCOZA is based on an approximate closure 
of the exact Ornstein-Zernike (OZ) equation that relates the two-body correlation function 
h(r) and the direct correlation function c(r): 



v{r) 



h(r) = c(r) + p Jd 3 r' c(|r - r'|) h(r') 



(6) 



where p = N/V is the number density of the fluid, V being the volume, and N the particle 
number. The SCOZA closure reads 



(7) 



h(r) = — 1 r < a , 

c{r) = cns{r) - K(p,f3)w(r) r > , 

where K(p, (3) is a state-dependent function of the density p and the inverse temperature 
f3 = 1/(&bT), and CHs( r ) is the direct correlation function of the hard-sp 



is assumed to be known, e.g. by the Verlet-Weis 



rere fluid. The latter 
19| or the Waisman 2(J parametrization. 



Here the Waisman parametrization has been used. Equation ([7]) is similar to the well-known 
optimized random phase approximation (ORPA) jjj: in particular, the requirement on h(r) 
for r < a, usually referred to as the core condition, is exact because of the hard core of 
the potential (j5J), while the expression of c(r) is clearly an approximation, since the off-core 
contribution to c(r) is assumed to depend linearly on the interaction w(r), hence having the 
same range as w(r) itself — the so-called OZ ansatz. However, in the ORPA the amplitude 
of w(r) is set to K = /3, while in the SCOZA K is a priori unknown, and must be determined 
so as to achieve consistency between the compressibility and internal energy routes to the 
thermodynamics. This requirement is expressed by the following condition on the isothermal 
reduced compressibility Xrod and the excess internal energy per particle U : 

d f 1 \ d\pU) 

~ P fl„9 ' ( 8 ) 



d(3 \XredJ dp 2 

where it is understood that Xrcd and U are determined respectively by the compressibility 
and internal energy routes, namely: 

1 



Xred 



1- pjd 3 rc{r) , (9) 
U = ~pjd 3 r[h(r) + l]w(r). (10) 

Equation would in fact be a trivial identity in a hypothetical exact description of the 
system, but this is not true anymore once an approximate closure such as Eq. ([7]) is intro- 
duced. Together with Eqs. JQ^M, (flUl). Eq. © yields a PDE for K(p,/3), that must be 
solved numerically. As in [l2|, Il8 |. Eq. (jSJ) is rewritten using the internal energy per unit 
volume u = pU as the unknown function: 

where we have set 

D(p,u) = ^(^-) . (12) 

The finite-difference algorithm used for the numerical integration of Eq. ffTT]) has been de- 
scribed in detail in Appendix [A] Here we observe that, irrespective of the specific dis- 
cretization of the PDE, a necessary step of the integration scheme consists in solving the OZ 
equation ([6]) supplemented with the closure (J7J). This amounts to solving the ORPA for an 
effective inverse temperature (3 e s = K. In several previous applications of the SCOZA, w(r) 



was chosen CIS db Yukawa tail |lllil4j]. a 
similar forms such as the Sogami-ESE 



inear combination of Yukawa tails 2l|, |22| , or other 
23] or related potentials 24J. If chs(^) for r > a is 



assumed to be vanishing as in the PY equation for the hard-sphere fluid [3[ , or is also given 
by a Yukawa tail as in the Waisman parametrization, these expressions of w(r) allow for an 
almost fully analytical solution of Eqs. (JHJ) , (El) - Although the integration of the PDE must 
always be carried out numerically, this leads nevertheless to a considerable simplification of 
the whole procedure. 

For the SW potential considered here, on the other hand, it is necessary to solve also 
Eqs. ([()]), (JZJ) in a fully numerical fashion. Moreover, as already found in 18| and recalled 
in the Introduction, for short-range wells the solution of the PDE is very sensitive to the 
position of the upper boundary p of the density interval, unless this is moved to very high 
values. Therefore, it is of paramount importance that the numerical algorithm employed 



remains convergent in this regime. In 



Eqs. ([6]), ([7]) were solved by the Labik, Mali- 
jevsky, and Vonka algorithm |25j]. However, in that work it was found that the algorithm 
failed to converge at high density and low temperature, so that for the shortest-range wells 
investigated there, namely for 0.25 < 5 < 0.5, the high-density boundary po was not moved 



beyond p a 3 = 1.15 
nally developed in 
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rerefore, here we have adopted a different method, which was origi- 



271 ] for the ORPA, and again applied to the ORPA for narrow SW 



and square-shoulder potentials in 



281 ] . This is based on the well-known property |3j that im- 



posing the core condition in the ORPA is equivalent to requiring the Helmholtz free energy 
^4rpa given by the random phase approximation (RPA) [3[ to be stationary with respect to 
variations inside the hard core of the quantity </>(r) defined as 

(f>(r) = c(r) - c H s(r) . (13) 

We recall that outside the core one has 0(r) = —(3w(r) in both the RPA and the ORPA. 

Specifically, the expression of A RP \ is 

(3A RPA (3A ns 
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where Ahs is the Helmholtz free energy of the hard-sphere fluid, and /hta and /Ri ng are 
given by 

/hta = ^yd 3 r0(r)[l + /i H s(r)] 
d 3 k 
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(15) 
(16) 



where hus( r ) an d <Shs(&) ar e respectively the two-body correlation function and the struc- 
ture factor of the hard-sphere fluid, and the hats denote Fourier transforms. In the equations 
above, /hta refers to to the high-temperature approximation, whereby the two-body cor- 
relation function h(r) is replaced by that of the hard-sphere fluid, while /rung refers to the 
sum of the ring diagrams in the expansion of Arpa in powers of 0(r) {3]. If fnmg is regarded 
as a functional of d>(r) inside the hard core, one finds 



Siting _ pA/l(r); 
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where we have set Ah(r) = h(r) — hns(r), and S(k) is the RPA structure factor given by 

S(k) = SUS(k) ; ■ (19) 

l- P S m (k)tf>(k) 

Note that, since hns(r) satisfies exactly the core condition, /hta does not depend on <p(r) 
for r < a, and therefore does not contribute to the functional derivatives of Arpa for r < a. 
Equation ffTTj) shows that the core condition is satisfied if and only if the RPA free energy is 
stationary with respect to variations of 0(r) inside the core, and Eq. (fIBI) shows that /pang 
is convex, so that this stationary point is unique and corresponds to a minimum of /Ring- 
As a consequence, one can solve the ORPA by minimizing /Ring numerically with respect 



to 0(r) in the core region [26j. Since Eqs. ( IT71) . ( |T8l) hold irrespective of the form of </>(r) 
outside the core, this procedure can be equally well applied to the present case where (see 
Eq. Q) 4>{r) for r > a is given by —K(p, (3)w(r). Therefore, in order to solve the OZ 
equation fl6]) with the closure (J7J), we can, for any given K, minimize the functional f^a ng of 
Eq. (fT6l) with <f>(r) = —Kw(r) for r > a. The difference with respect to the ORPA consists 
in the fact that in the ORPA, the minimum of the functional gives the Helmholtz free energy 
of the system as predicted by the energy route, while this is not true anymore for a generic 
K. In fact, one finds 

= -pj dh Ah(r)w(r) = -2{U- U U ta) , (20) 

where U is given by Eq. ( fTUl) . and Uuta is the corresponding quantity in the high-temperature 
approximation: 

Ukta = \pj d3 r [^Hs(r) + %(r) . (21) 



For K = (3, Eq. ( 120]) coincides with the usual thermodynamic relation d(/3A/N)/d/3 = 
U with A given by Eq. ( 1T4"|) . so that A is indeed the Helmholtz free energy obtained by 
integrating the internal energy U with respect to (3. On the other hand, this is not the case 
if K has a nontrivial dependence on (3. Therefore, unlike in the ORPA and the RPA, in the 
SCOZA Eqs. (HU), (US}, (USD do not give the Helmholtz free energy Nevertheless, Eq. flU} 
turns out to be very useful to overcome a difficulty related to the numerical solution of the 
SCOZA. Specifically, the calculation of the quantity D(p,u) defined in Eq. ffT2~]) requires 
that Eqs. fl6]), ([7]) must be solved at fixed p and U, rather than at fixed p and K, i.e., for 
any given density one needs to find the value of K such that the internal energy per particle 
takes a given value U = U. One way to do this would be to change K by trial and error, 
solve Eqs. ([H]), (EJ) with respect to <f)(r) inside the core for each K, and find the corresponding 
U by Eq. ([TO]) , until the condition U = U is met up to a prescribed accuracy. However, we 
found that one can solve Eqs. (JBJ), (J7|) simultaneously with respect to <p{r) and K in a single 
optimization run. To this end, let us observe that one has 

<9 2 /Ring p f d 3 k 



-pjd 3 vAh(r)w(r) + AU, (24) 



dK 2 2 J (2ttY 

where S(k) is given by Eq. f lT9|) . We now introduce the Legendre transform of fm ng 

S = /R in g - = /Ring + 2K(U - C/hta) • (23) 

In the case of the ORPA in which K = (3, S coincides with twice the excess entropy per 
particle with respect to the hard-sphere fluid. Equations ( )20|) . (122]) imply that, if 2(U— Uuta) 
is fixed at some given value AU, S is a convex function of K such that 

dS 
dK 

so that its minimum is obtained for that K which, when inserted into Eq. (J7]), will yield AU 
via Eqs. (TiTJ]) . (12~T]) for the internal energy. Therefore, in order to solve the OZ equation ([6]) 
with the closure (J7]) for a fixed density p and internal energy U, it is sufficient to minimize 
S, regarded as a function of K and a functional of 0(r) for r < a. The numerical procedure 
adopted for the minimization is described in detail in Appendix [B] 

A difficulty related to the aforementioned request that Eqs. (EJ), (EJ) must be solved also 
at very high density, is that in this regime the structure factor of the fluid develops very high 
and narrow, "pseudo-Bragg" peaks. These do not have a direct physical meaning, since they 
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occur in a density range where the stable phase would actually be the solid one. Nevertheless, 
they must be taken into account when solving Eqs. ([6]), ([7]) numerically. In fact, this problem 
appears already for the simple hard-sphere fluid: if one starts from the (analytical) chs(^) 
given by the Waisman parametrization, finds hns(k) by the OZ equation and obtains 
^Hs( r ) by performing directly the inverse Fourier transform of hns(k), it is found that for 
p* > 1.2, the core condition is very poorly satisfied. The reason for this is that many 
of the peaks of /ihs(&) will simply be missed by the discrete sampling of the wave vector, 
unless the number of point used in the transform is so large (about 2 18 ), that the numerical 
computation becomes unpractical. In order to cope with this problem, the contribution of 
these peaks to the inverse Fourier transform has been determined analytically by the method 
described in Appendix O 

The initial condition required by the PDE (jTTJ) at /? = and the boundary conditions at 
p = and on the spinodal curve, i.e., the locus of diverging compressibility below the critical 
temperature, have been detailed elsewhere 12J. The high-density boundary condition is more 
problematic in this context, since this is not known exactly while, as stated in Sec. [H the 
results for narrow SW are quite sensitive to both the position of the high-density boundary 
Po, and the kind of approximation used for u at p = po- This point will be considered in 
more detail in Sec. IIII( here we remark that in this investigation the high-density boundary 
condition has not been imposed on the second derivative of u with respect to the density 



d 2 u/dp 2 as in a number of previous works 12H14I. [18} , but rather directly on u, because we 



found that the results obtained in this way were less sensitive to the choice of w(po, (3). 

The vapor-liquid coexistence curve is obtained by imposing the conditions of thermo- 
dynamic equilibrium, i.e., by finding the densities p v , pi which give the same pressure P 
and the same chemical potential p on the liquid and vapor branches of the sub-critical 
isotherm. P and p are determined by integrating with respect to (3 the exact identities 
d(/3P)/d/3 = — u + pdu/dp, <9(/3p)/<9/3 = du/dp, which are fulfilled also by the SCOZA 
because of the consistency condition (jSJ). 

Concerning the parameters used in the numerical calculation, we typically set the density 
step at Ap* = 10 -3 , and the initial inverse-temperature step at A/3q = 10~ 3 . In order to 
locate the critical point accurately, A/3 is decreased as the critical point is approached, and 
then gradually expanded back to A/?o below the critical temperature. We found that the 
results are not very sensitive to the density or temperature mesh, e.g., even a substantial 
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coarsening of the mesh such as Ap* = 10~ 2 , A/3q = 10~ 2 leaves the results nearly unaffected, 
except for the coexistence curve in the neighborhood of the critical point. In this case, 
a very high accuracy in the numerical calculation is required to fulfill the conditions of 
thermodynamic equilibrium, so that a coarse mesh results in a large gap at the top of 
the coexistence curve. On the other hand, the numerical output turns out to be much more 
sensitive to the discretization used in the numerical Fourier transform. We then used a rather 
small step in real space, Ar* = 5 x 10 -4 , and a large number of mesh points, N = 2 15 . We 
employed the fast Fourier transform routines developed by Takuya Ooura at the Research 
Institute of Mathematical Sciences, Kyoto University, and made freely available online. 

Whenever the function to be transformed is discontinuous, such as Ah(r) or Ac(r), 
it is best to subtract off the discontinuity, and add the slowly decaying transform of the 
discontinuous part determined analytically at the end. Here we implemented this procedure 
on both the function and its derivative, so that the numerical transform is done on a smooth 
function that is not only continuous, but also differentiable. If Df(ro), Df'(ro) denote the 
discontinuities respectively of / and its derivative /' at r = r , i.e., Df(r Q ) = \im r _^ r - f{r) — 
lim r ^ r + /(r), Df'(r ) = lim r _ >r .- fir) — lim r ^ r + f'(r), the function / that is transformed 
numerically is given by 

/>) = f(r) - [Df(r ) + (r - r )Df'(r )} [1 - 6(r - r )] , (25) 

where 0(r) is the Heaviside step function defined by 0(r) = 1 for r > 0, O(r) = for r < 0. 
For the SW potential studied here, the discontinuities of Ah(r) and Ac(r) are located at 
r = a and r = er(l + 5). 

The accuracy within which the core condition is fulfilled can be determined a posteriori 
by calculating the radial distribution function g(r) = h(r) + l inside the repulsive core r < a. 
Typically, we found that for T ~ T c and high densities such as p ~ 0.9er~ 3 , the order of 
magnitude of g(r) is g(r) ~ 10~ 8 -10~ 7 for O.lcr < r < a and g(r) ~ 10 _6 -10 -3 for r < 0.1a. 



III. RESULTS 

We have employed the numerical solution of the SCOZA described in Sec. [Ill to study 
SW potentials of several widths 5. For large 5, our results are basically identical to those 



obtained in 



181 ] . As can be inferred from the comparison between Tab. [J shown further in 
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this Section and Tab. I of 



18| . the critical points for S > 0.5 differ at most by ~ 0.2%. Here 



In this 



we discuss in detail our data for 5 < 0.25, which is the smallest value studied in 
regime, the liquid-vapor transition considered here is actually expected to be metastable 
with respect to freezing 
As already found in 



and anticipated in Sees. HJ [TTJ we also have observed a strong 
influence on the results of the high-density boundary condition, that becomes more and more 
important as S decreases. This is especially true for the critical density p c and the high- 
density branches of the spinodal and coexistence curves. In order to assess this sensitivity, 
the high-density boundary p was initially set at pj$ = 1, and then moved forward to values 
as high as p* = 1.4, which is very close to the close-packing density p* Q = \/2. Moreover, 
for each po several different approaches for the boundary condition were tried. Specifically, 



the interna 
ORPA" 



energy w(p ,/3) has been determined by the HTA; the ORPA; the "non-linear 



161 ]. whereby /3 is replaced by the Mayer function e?* — 1 as the amplitude of the 



potential in c(r), i.e., one sets K* = — 1 in Eq. ([7]); the EXP [3[, which is obtained by 
setting gEXp(r) = gns( r ) exp[A/?,oRPA (?")], where g{r) is the radial distribution function, and 
A/iorpa(^) = ^orpa(?") _ ^hs(t); the linearized version of EXP known as LIN [3[, where one 
has guN( r ) = fi'Hs( r )[l + A^orpa (?")]; an d, f° r ver y narrow wells, the solution of the PY 
equation in the SHS limit Q . For the values of 5 discussed below, the results were found to 
be enough stable with respect to a further increase of po or a change of the approximation 
used for w(po,/3) to warrant comparison with other predictions. 

In Fig. [1] we have reported the SCOZA coexistence and spinodal curves for 5 = 0.25 
obtained with Pq = 1.2 and u(po, 0) given by the non-linear ORPA. The critical temperature 
and density predicted by the SCOZA are T* = 0.769, p* = 0.355. These values are close, but 



not identical to those found in [18| for the same value of 5, namely T* = 0.761, p* = 0.343. 
The difference is to be traced back to both the somewhat lower value of po used there, 
Po = 1.15, and the fact that the boundary condition at po was imposed on d 2 u/dp 2 rather 
than u. As observed in Sec. HIl we found that this enhanced the sensitivity of the results 
to the choice of the approximation used at po, namely, the HTA in the case of 18|. The 
main difference between the two coexistence curves is that the liquid branch of the present 
one is slightly shifted to higher densities with respect to that obtained in [l^ . Nevertheless, 
the two curves are close to each other, and the comparison with the simulation results 



29 



34j shown in Fig. [T] confirms the trend already found there, i.e., the liquid branch of the 
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FIG. 1: Liquid-vapor coexistence curve and critical point of the SW fluid of width 5 = 0.25 in the 
density-temperature plane. Solid line: SCOZA. Full circle: SCOZA critical point. Squares with 
error bars: simulation data by Vega et al 



Tri ang les: simulation data by Liu et al. 



301 ] . Diamonds: simulation data by Elliott et al. 31] • 



al. 



291 ] Inverted triangles: simulation data by del Rio at 



32|. Circles: simulation data by Orea et al. 



34i | . Crosses: simulation data by Pagan et al. 



33(]. 



The SCOZA spinodal has also been shown by the dashed line. 

SCOZA coexistence curve underestimates the simulation results. The different simulations 
reported in the figure are enough consistent with one another to indicate that this is indeed 
an inaccuracy of t he SCOZA. The critical temperature is in good agreement with most 
simulation results 



29 



30 



32H34j . while the critical density is underestimated with respect 
to all of them, although it lies within the error-bar given in jsoj] , and the discrepancy is 
generally smaller than that observed on the liquid branch of the coexistence curve. 

Figures [21 [3] show the SCOZA spinodal and coexistence curves respectively for 5 = 0.15 



and 5 = 0.1, together with simulation data for the coexistence curve and the critical point 
29 



10 



33 



361 ] . The SCOZA results were again obtained using the non-linear ORPA as the 



boundary condition at po, and setting p$ = 1.4. Unlike what found for 5 = 0.25, the 
coexistence curve of the SCOZA is wider than most simulation predictions. This behavior 
is largely due to the fact that now the SCOZA critical temperature is appreciably higher 
than the simulation ones. 
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FIG. 2: Same as Fig.[[]for 5 = 0.15. Solid line: SCOZA. Full circle: SCOZA critica 



simulation data by Orea et al 
simulation data by Liu et al 



34i | . Crosses: simulation data by Pagan et al. 



point. Circles: 



33J. Triangles: 



29(]. Dashed line: SCOZA spinodal. 



As a general remark, the overall agreement with the phase diagram predicted by the 
SCOZA appears to be worse than that found for HCY potentials of comparable range 
151 ]. although for 5 = 0.15 and 5 = 0.1 a quantitative assessment of the performance of the 
SCOZA is difficult, because of the rather large discrepancies shown by different simulations. 
A perhaps more conclusive comparison can be made with the recent simulation results for 
the critical point obtained in 10|], some of which have been reported in Tab. [I] together with 
the SCOZA predictions. These indicate that the SCOZA for the SW fluid underestimates 
p c , the discrepancy on p c being larger than that on T c . In a previous investigation of the 



SCOZA for narrow HCY fluids 15J, p c was again found to deviate from simulations more 
than T c , although in that case p c was overestimated by the SCOZA. 

A qualitative feature common to the SCOZA phase diagrams shown in Figs. [H, EJ |3] is 
that the coexistence curve is much wider then the spinodal curve, so that even at densities 
near the critical one, and a fortiori at low or high densities, there is a sizable temperature 
interval where the fluid may survive in a single phase as a metastable state, before spin- 
odal decomposition sets in. According to this, in the experimental practice one should be 
circumspect before using spinodal decomposition to identify equilibrium phase separation 



15 



H 0.45 




0.2 0.4 0.6 0.8 



FIG. 3: Same as Fig. Q] for 5 = 0.1. Solid line: SCOZA. Full circle: SCOZA critical point. Squares 
with error bars: simulation data by Duda 361. Circle with error bars: simulation critical point by 

n □ . 

Skibinsky et al. [351 ] . Diamond: simulation critical point by Largo et al. [10(. Dashed line: SCOZA 



spinodal. 



on the assumption that the coexistence and spinodal curves run close to each other. In 
fact, according to the SCOZA they do not, especially at the low volume fractions that are 
frequently considered in experiments on colloidal dispersions. 

For wells such that 5 < 0.05, the dependence of the SCOZA on the high-density boundary 
condition becomes so strong, that the results do not show any clear trend towards saturation 
in the region < 1.4. It is possible that this could be achieved by moving p to even higher 
values, bigger than that corresponding to the physical close packing. However, with the fully 
numerical algorithm considered here we failed to obtain solutions of Eqs. ([6]), (IT)) beyond 
close packing, because in this regime the correlation functions become very singular, and 
even the prescription described in Appendix lAl proved insufficient, at least for the number 
of points in the numerical Fourier transform that can be handled in practice. Therefore, the 
results that we have obtained for 5 < 0.05 cannot be considered meaningful. 

It is worthwhile pointing out that the sensitivity on the boundary condition at high 
density is much stronger for the SW potential considered here than for the HCY potential. 
In order to compare the two cases, it is useful to resort to the NF mapping mentioned in 
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SCOZA MC 



U 


± c 


He 


M c 


n* 


0.8 


1.953 


0.247 


1.940 


0.263 


0.5 


1.211 


0.272 


1.220 


0.310 


0.3 


0.850 


0.328 


0.847 


0.376 


0.2 


0.683 


0.382 


0.667 


0.421 


0.1 


0.493 


0.459 


0.4780 


0.478 



TABLE I: Critical temperature and density of the SW fluid for several values of the well width 
5, as obtained by the SCOZA and Monte Carlo (MC) numerical simulations performed in [lo| . 



According to the estimate given in 



101 ] for the case 5 = 0.05 (not reported here), the errors in the 



MC critical parameters are about ±0.14% for T* and ±1.56% for p*. 

Sec. HI according to which a fluid of hard-core particles of diameter a with a given attractive 
tail potential at reduced temperature T* can be mapped into a SW fluid with the same 
a, and a range determined so that its second virial coefficient B 2 {T) coincides with that 
of the original fluid for T = T*. We remark that here this mapping is used without any 
assumption as to the behavior of the system for interactions of vanishing range, so its use is 
fully legitimate also within the SCOZA. A HCY fluid with inverse-range parameter z = 5.5 
at T* = 0.48, which is the critical value according to the SCOZA, should then be equivalent 
to a SW fluid with 5 = 0.1 at the same T*. The SCOZA critical temperature of this SW, 
T* = 0.49, is indeed close to that of the HCY, in agreement with the prediction from the 
mapping, while the difference between the critical densities, p* = 0.459 for the SW and 
p* = 0.427 for the HCY, although larger, is still below 10%. However, for the HCY fluid 
setting Pq = 1 is sufficient to obtain a stable critical point, while for the SW must be set 
at a substantially higher value, Pq = 1-4- 

Further insight into the difference between the two potentials is obtained by checking the 
behavior of the effective range 5 of the HCY below T c : we find that, as the temperature 
decreases, the effective range also decreases. As a consequence, for T < T c it should be 
possible to map the densities at coexistence of the HCY fluid into those of a SW fluid at 
the same T*, but with a lower T*, i.e., at a higher rescaled temperature T/T c . This implies 
that, according to the NF mapping, the coexistence curve of the HCY fluid in the T-p 
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FIG. 4: Solid line: SCOZA coexistence (outer) and spinodal (inner) curves for a SW fluid of width 
5 = 0.1. Dashed line: same for a HCY fluid with inverse-range parameter z = 5.5. The upper panel 
shows the curves in the temperature-density plane, with both variables rescaled by their critical 
values. In the lower panel the rescaled temperature has been replaced by the rescaled stickiness 
parameter r/r c Q|. Note how the curves for the two potentials get closer in the latter case. 

plane should be narrower than that of the "equivalent" SW at T*. That this is actually 
the case is shown in the upper panel of Fig. HJ where the coexistence and spinodal curves 
of the HCY fluid with z = 5.5 and the SW fluid with 5 = 0.1 have been plotted using 
rescaled temperatures T/T c and densities p/p c - It is evident that the SW coexistence curve 
is appreciably wider, especially when the liquid branch is considered. It is also worthwhile 
considering the behavior of the coexistence and spinodal curves when the temperature is 
replaced by B 2 (T), or equivalently by the stickiness parameter r obtained by replacing I?f HS 
with B 2 (T) in Eq. ([3]). The NF mapping predicts that, with this choice of the temperature- 
like variable, the phase diagram of narrow attractive potentials should become universal. 
The two fluids considered here do give similar values for r at the critical point, namely, 
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FIG. 5: Reduced second virial coefficient of the SW fluid at the critical temperature B^Tc) as a 
function of the well width 5. Squares: SCOZA. Circles: simulation results by Largo et al. [lo| . 
The dotted line is a guide to the eye. 

r c = 0.115 for the SW and r c = 0.107 for the HCY. If the rescaled temperature T/T c is 
replaced by the rescaled stickiness t/t c , as in the lower panel of Fig. HJ the vapor branches of 
the coexistence and spinodal curves become nearly coincident, and the discrepancy between 
the liquid branches of the coexistence curve becomes substantially smaller, although the 
SW still gives a wider curve. Quite curiously, the situation with the liquid branches of the 
spinodals is reversed with respect to that shown in the upper panel, 1.6., clS cl function of r 
the SW spinodal is narrower than the HCY spinodal, while as a function of T it is wider. On 
the basis of the above considerations, the sensitivity of the SW potential to the high-density 
boundary condition may depend at least in part on the fact that its phase diagram is shifted 
towards higher densities than those of a HCY fluid of equivalent range. 

The behavior of B 2 (T C ), shown in Fig. [5j is in semi quantitative agreement with the 



simulation results 10] . In both cases, the B 2 (T C ) increases as 5 decreases. However, because 
of the impossibility to obtain meaningful results for the critical point of narrow wells such 
that 5 < 0.05, we are not able to say what the fate of the SCOZA B 2 (T C ) will be in the 

nn 

limit 5—7-0, namely, if it will tend to a finite value as predicted by the simulations p, I10| . 
The arguments put forth in Sec. HV] rather suggest that the SCOZA would not yield a finite 
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FIG. 6: Compressibility factor /3P/p (upper panel) and excess internal energy per particle U* 
(lower panel) of the SW fluid with 5 = 0.1. Solid line: SCOZA for T* = 1. Dashed line: SCOZA 
for T* = 2. Circles: simulation results at the same temperatures \s7\. 

value of E>2(T C ) for vanishing attraction range. 

We now briefly consider some thermodynamic properties of the SW fluid. Figure E] 
displays the compressibility factor flP/p and the internal energy per particle U* for 5 = 0.1 
along the isotherms T* = 2 and T* = 1, compared with simulation results [3?]]. The 
agreement is very satisfactory both at low and high density, although it is fair to point out 
that the temperatures at which the simulations were performed can be considered quite high, 
since even T* = 1 is still more than twice the expected critical value reported in Tab. [B A 
comparison encompassing temperatures closer to T c is shown in Fig. [3, where (3P/ p and U* 
for 5 = 0.25 are shown along the isochore p* = 0.42 as a function of the inverse temperature 



(3, and compared with simulation results 38[. The agreement is again quite good down 



to the lowest temperature attainable by the SCOZA, i.e., that at which the isochore hits 
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FIG. 7: Compressibility factor j3P/ p (upper panel) and excess internal energy per particle U* 
(lower panel) of the SW fluid with 5 = 0.25 and p* = 0.42. Solid line: SCOZA. Circles: simulation 
results 



the high-density branch of the spinodal curve. Inspection of Fig. 2 of 38J suggests that 
the simulation data at the lowest temperatures, roughly for j3 > 1.33, would correspond to 
unstable states in the thermodynamic limit. 

It is also worthwhile examining how the present theory performs for the correlations 
of the fluid. Figure [S] shows the structure factor S(k) for a SW with 5 = 0.1 and two 
different states, one at moderate density p* = 0.5 and relatively high temperature T* = 2, 
and another at lower density p* = 0.229 and temperature T* = 0.5 near the critical value. 

39l. The 



In both cases, the SCOZA reproduces quite accurately the simulation results 
correlations in real space are considered in Fig. |9j where the radial distribution function 
g(r) predicted by the SCOZA, again for 5 = 0.1 and T* = 0.5, is compared with simulation 



data 



40j at the densities p* = 0.2, p* = 0.4, p* = 0.8. At low density, the theory agrees 
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FIG. 8: Structure factor of the SW fluid with 5 = 0.1 for T* = 2, p* = 0.5 (upper panel) and 
T* = 0.5, p* = 0.229 (lower panel). Solid line: SCOZA. Circles: simulation data by Shukla et 



al. 



39l | extracted from their Figs. 6 and 8. 



closely with the simulation, both inside and outside the well region. This is not surprising 
since, as pointed out in Sec. [U the SCOZA for the SW potential becomes exact at low 
density, while this is not the case for other tail interactions. On the other hand, as p is 
increased the SCOZA overestimates both the contact value g{cr + ) of g(r), and the amplitude 
of the discontinuity at the well boundary £ = (1 + 5)cr. This also testifies the quite special 
status of the SW potential within the SCOZA: in fact, for the HCY potential 13j the 
SCOZA was found to underestimate g{<J + ), as one would expect from a theory whose c(r) 
is linear in the perturbation potential. In the present case, the behavior at high density is 
instead reminiscent of that of the "non-linear ORPA" studied in 16|. Another feature worth 



attention is the hump displayed by the simulation g(r) for p* = 0.8 and r* ~ 1.8, which is 
the precursor of the discontinuity for r* = a/3 displayed by g{r) in the sticky limit a a 
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FIG. 9: Radial distribution function of the SW fluid with 6 = 0.1 for T* = 0.5 and p* = 0.2 (top 
panel), p* = 0.4 (central panel), p* = 0.8 (bottom panel). Solid line: SCOZA. Circles: simulation 



data by Largo et al. 



40l | extracted from their Fig. 6. 



This is not reproduced by the SCOZA, most likely because of the assumption implied by 
Eq. ([7]) that c(r) has the same range as the potential. 

A closer look at the g(r) at contact and at the well boundary is given in Fig. [T0J where 
g(<J + ), g{£,~), g{£ + ) for 5 = 0.1 are shown along the isotherms T* = 0.5 and T* = 1. As one 
may expect, the aforementioned discrepancy with the simulation data that sets out at large 
p becomes less severe as T is increased, the results for T* = 1 being in good agreement with 
the simulation 40] even at the highest p considered. 
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FIG. 10: Radial distribution function of the SW fluid with 5 = 0.1 for r — > cr + , r — > r — > £ + , 
where £ = (1 + 5) a is the value of r beyond which the potential vanishes. Squares, triangles, 
and circles represent the SCOZA results for g(a + ), g(£,~), and <?(£ + ) respectively. Crosses, pluses, 



and starbursts represent the simulation data obtained by Largo et al. 



40| for the same quantities. 



Dotted lines are a guide to the eye. Upper panel: T* 
isotherm. 



1 isotherm. Lower panel: T* = 0.5 



IV. THE STICKY LIMIT OF THE SCOZA 



As discussed in Sec. IIII[ the numerical solution of the SCOZA for the SW potential 
developed here becomes unpredictive for wells such that 8 < 0.05, still considerably wider 
than the HCY potentials that were considered in [l^, Tsj] . However, in the limit <5 — it is 
possible to study the behavior of the SCOZA by adopting a different standpoint, namely, by 
imposing thermodynamic consistency on the solution to Eqs. ()6]), fl7J) after the limit 5^-0 
has been taken. The clue to this procedure lies in the observation made in 42| that in the 
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limit 5 — > 0, any OZ closure like Eq. (J7J), such that one has /i(r) = — 1 for r < a and c(r) = 
for r > (1 + <5)<t, gives a fixed form for c(r): 

c (r) = |_^- a 2 + r/ [6( a + 6) 2 - aX] r* - ^r* 3 j 9(1 - r*) + ^ 5(r* - 1) , (26) 

where 0(x) and 5(x) are respectively the Heaviside and Dirac functions of argument x, and 
we have set rj = irpa 3 /6. The quantities a and b are given by 

1 + 2r] \r] 



[I-t]) 2 1-v' 
3rj Xt] 



(27) 
(28) 



2(1 - V y 2(1-77)' 

Therefore, in the limit 8 — > all the specific features of the OZ closure considered are 
conveyed into a single, state-dependent parameter A. The latter is related to the stickiness 
parameter r defined in Eq. flH]) and to the contact value of the cavity function y(r) = 
g(r) exp [/3v(r)] by the expression 

\=*W. (29) 
r 

By applying the compressibility rule (j9]) to Eq. ( 1261) . one obtains 

1 



Xred 



a 2 , (30) 



whereas the energy route gives 43] 



pa_p^_ r + - dT ,M^i (31) 



iV N ' J T t' 

If the form of A(r/, r) obtained by the PY closure is substituted into Eqs. (l2"6]) -f l3~T|) . the 

aforementioned solution of the SHS model due to Baxter is recovered. On the other 
hand, as discussed in detail in 



42], different choices for A are also possible. Since one 
of the most problematic features of the original SHS model is its lack of thermodynamic 
consistency, it is tempting to determine A via the SCOZA, thus enforcing consistency between 
the compressibility and energy routes given by Eqs. fl30l . ( 13~T|) . If we introduce the quantities 

A = rfX , (32) 
t = - , (33) 

T 



the SCOZA equation reads 



^ x c> A d 2 A 

^,t)t- = -, (34) 
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where V(rj,t) is given by 



V(r],t) 



77 (1 + 277) - (1 -77)A 



(35) 



The PDE ([34") c an be solved numerically by the finite-difference scheme described in the 



Appendix of 



44J, or the similar one given here by Eq. flAlj) ; unlike in the case of non 



vanishing S, the diffusion coefficient is given explicitly by Eq. ( )35l) . so there is no need to 
evaluate it by finite differences as in Eq. ( 1A2I) . As usual, for the problem to be completely 
specified, Eq. ( I34"l) must be complemented with the initial and boundary conditions. As for 
the initial condition, we observe that the numerical integration cannot be started exactly 
from t = 0, because for this value of t the diffusion coefficient 1/(£D) would become infinite. 
Instead, the integration was started from a small, but non null t > 0, using as the initial 
condition for A the leading term in an expansion in powers of to 



Mv,to) - t ip(ri) . 



(36) 



By substituting Eq. f l36|) into Eq. fl34"|) we find that the //-dependent term 99(77) must satisfy 
the ordinary differential equation 1 



rfV_ 2 1 + 2// 

drf 77 2 (1 — 77) 3 

This is solved numerically with the conditions 



d 2 cp 







r)=0 



dr] 2 



(37) 

(38) 
(39) 



11=0 



that follow from Eqs. (129]) , (l3"2l. (j35jl and the property |42j 



limy((x) 



1 . 



(40) 



We have verified a posteriori that the numerical solution of Eq. (!34|) is virtually independent 
of the particular choice of to, as long as the latter is sufficiently small. Equations (129]) . (132|) . 
also imply for the low-density boundary condition 



A( v = 0,t) = 



for every £. 



(41) 



1 Although for r — > 00 the contact value of the radial distribution function coincides with the PY result, 
y(a) docs not because of the double limiting procedure 5 —> 0, r — >• 00 in the prefactor e' 9 "^ cr+ '. 
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As observed in Sec. [ITJ, the choice of the high-density boundary condition entails a certain 
amount of arbitrariness. A natural possibility is using for A the expression Aeax obtained by 
the analytical solution of the SHS model within the PY approximation j^j]: 

A(?7o, t) = J7oA Ba x(?7o, t) for every t, (42) 

where rj denotes the position of the high-density boundary. Quite disappointingly, although 
perhaps not surprisingly, the outcome of the numerical integration of Eq. (I3"4"|) with the 
initial and boundary conditions given by Eqs. ( 13T)]) . (141 p . fT4"2"|) suffers from the same problem 
experienced in Sec. IIHI for wells of non vanishing width, namely, a strong dependence of 
the results on the position of 770- This is illustrated in Fig. [TTJ where the inverse reduced 
compressibility 1 / Xrcd obtained for different choices of 770 has been plotted along the isotherm 
r = 0.11, which is close to the critical one according to numerical simulations For the 
solution to be deemed reliable, the various curves, at least for 770 greater than a certain 
threshold, ought to lie near one another — a requirement that is manifestly not satisfied. 
Some of the isotherms are not even consistent with r = 0.11 being close to the expected 
critical value. We have also verified that this state of affairs is not changed by different 
choices of the high-density boundary condition. 

We can, however, resort to a different method to solve Eq. (134"]) . that does not require 
the high-density boundary to be specified. This is made possible by the fact that, unlike in 
the case of wells of non vanishing width discussed in Sees. [TTl IHIl the diffusion coefficient of 
Eq. fl34"l) is known explicitly as a function of r\ and r. Let us assume that A can be expressed 
as a power series in t 

00 

A(r ] ,t) = Y,t n ou n (r ] ), (43) 

n=l 

where the zeroth-order term has been omitted, since A must vanish for t — > 0. By inserting 
Eq. (J43J) into Eq. (134|) and equating the terms of the same degree in t, Eq. (134)) is replaced 
by an infinite set of ordinary differential equations: 

d 2 u n l + 2r/ 2 ^. 

which gives back Eq. (1371) for n — 1. Since this set contains only backward references, it can 
be integrated numerically to any desired order n. Only the low-density behavior of u n is 
required. Specifically, one needs the values of u n , du n /dr], and d 2 u n /dr] 2 at r\ = for n < n. 
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FIG. 11: Inverse reduced compressibility 1/xred of the SHS fluid as a function of the reduced density 
p* as found via the SCOZA equation (|34p with a high-density boundary condition computed in 
the PY approximation 2fl. The different lines refer to the boundary packing fractions rjo = irp^/G 
given in the legend. The stickiness parameter is set at r = 0.11, which is close to the critical value 
according to numerical simulations Q]. 

If it is assumed that u n can be expanded in power series of rj, it is found that for n > 6 
all the above quantities are in fact vanishing, while for n < 6 they can be obtained from 
exact diagrammatic expansions available in the literature 
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45| . Unfortunately, the radius 



of convergence £r of the series (143)) corresponding to most values of the packing fraction rj 
comes out to be too small for it to be of direct usefulness. This is apparent from Fig. [T21 
which shows how £r. rapidly decreases as the density is increased. In particular, at the critical 
value of t predicted by simulations p* should be well below 0.2 to obtain convergence. 
This is again disappointing, yet not uncommon: a power series obtained as a formal solution 
of a differential equation may easily turn out to be divergent. In those cases, one can try a 
resummation technique in order to provide an analytic meaning to such a divergent series. 



To this end, we have employed Borel summation method |46j, whereby one starts from the 
exact equality 

00 -in , oo ^.+00 n 

AMHE^-i^HtE / dse-^-uM (45) 



28 



35 
30 
25 
20 

Qi 
■(— > 

15 
10 
5 


FIG. 12: Radius of convergence tp, of the power series (|43p . estimated by means of the ratio 
convergence test, as a function of the reduced density p* . The horizontal dashed line marks 
the value of the inverse stickiness t corresponding to the critical point according to numerical 
simulations Q]. 

and formally interchanges the series and the integral operators to obtain 

i ;-+oo 

AbM) = -/ dse-'ft&fas), (46) 
£ Jo 

where A(?7, s) is the regularized power series 

°° s n 

Hv,s) = J2 — UJ n(v), (47) 

71. 

n=l 

that has an improved radius of convergence because of the factorial in the denominator. 
Where the original series is convergent, the two functions A and Ab are in fact identical. 
The latter, however, yields a finite value in a much larger domain, providing the analytic 
continuation of the function defined by the original series. If one is willing to trust the 
Borel sum Ab in this enlarged region, then a puzzling result is found: the liquid-vapor phase 
transition disappears, at least for reasonable values of the packing fraction and temperature. 
In Fig. [13] the inverse reduced compressibility of the fluid 1 / Xrcd is shown as a function of 
t along the 7] = 0.25 isochore. While the compressibility route of the PY approximation 
predicts the compressibility to diverge around r ~ 0.08, the solution of Eq. fl3~4"|) evaluated 
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FIG. 13: Inverse reduced compressibility 1/xrcd of the SHS fluid as a function of the inverse 
stickiness t along the ij = 0.25 isochore. The compressibility route of the PY approximation 
(dashed line) predicts the compressibility to diverge at t ~ 12. According to the solution of 
Eq. (|34j) evaluated by means of the Borel summation method (solid line), the isochore does not 
cross the spinodal at all. 

by means of Borel summation method gives instead a finite Xred even if the temperature 
is lowered further. This result is qualitatively similar to that found by direct numerical 
integration of Eq. (I34j) when the high-density boundary rjo is pushed to (unphysically) high 
values, but is clearly at odds not only with the PY solution of the SHS model but also 



with the simulation results 



0. 



The present investigation supports the conclusion that in the 
limit of vanishing well width S — > 0, the SCOZA critical temperature decreases too rapidly 
to give a finite value of B 2 (T C ), so that B 2 (T) or, equivalently, the stickiness r, do not repre- 
sent a convenient way to measure the temperature in this limit. Imposing thermodynamic 
consistency on the SHS model does not lead to a better quantitative agreement with simu- 
lation compared to the highly inconsistent PY solution, but merely makes the liquid-vapor 
transition disappear. 
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V. CONCLUSIONS 



We have performed an investigation of the SCOZA for SW potentials, focused on the 
regime of narrow wells which is the most relevant for the modelization of colloidal particles. 
This study was prompted by the lack of a theory able to account quantitatively for the 
liquid-vapor phase transition in a system of sticky hard spheres: in the well known and 
widespreadily used PY approximation [2|, the various routes to the thermodynamics of the 
fluid are inconsistent with one another, and agree only qualitatively with the results of 
computer simulation hjl. On the other hand, for HCY potentials the SCOZA can be solved 

nn 

semi-analytically, and has already been usefully applied to narrow-range tails UAUM. The 
rationale for turning to the SW potential within the SCOZA is that this interaction has the 
unique feature of having the SCOZA become exact at low density, thereby recovering the 
correct B 2 (T) which plays such an important role for narrow interactions. 

However, the results presented here show that this peculiarity turns more into a liability 
than an asset. For wells of non vanishing width S, after extending the study already made 



in 



18] to (slightly) smaller 5, we cannot but confirm the conclusions of that work, namely: 



i) the results are found to be very sensitive to the boundary condition at high density that 
is needed to integrate the SCOZA PDE, actually much more so than in the case of HCY 
potentials of comparable range. Such an undesired feature can be tamed by pushing the 
density of the boundary to very high values, but this has the effect of making the numerical 
solution of the PDE unfeasible for 5 < 0.05. ii) In the range of 5 for which meaningful 
results can be extracted, the overall agreement with simulation data, at least for the phase 



13r|l5|, although 



diagram, is actually worse than what was found for the HCY potential 
a sharp assessment of the performance of the SCOZA is made difficult by the considerable 
discrepancies between different simulations. 

If the sticky limit is taken from the outset, and subsequently thermodynamic consistency 
is imposed, the SCOZA PDE is obtained explicitly in closed form, but the high-density 
boundary problem is again found to be looming. If the PDE is turned into an infinite 
set of ordinary differential equations and the low-density behavior alone is imposed, the 
liquid-vapor phase boundary disappears altogether. This supports the conclusion that as 
the attractive potential range vanishes, the SCOZA critical temperature will decrease more 
rapidly than what is necessary to yield a finite B 2 (T C ). 
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Stating that the SCOZA is inadequate for SW interactions in the sticky limit is merely 
an observation, and does not shed any light on the reason behind the fact. In this respect, 
one should recall the result of Stell js] mentioned in Sec. HI who showed that the SHS model 
is strictly speaking ill-defined if a class of problematic configurations of particles is not ruled 
out, for instance by introducing a small amount of polidispersity into the system. Can 
it be that the (possibly more accurate) SCOZA approach exposes the malicious effect of 
the diverging clusters, which are simply neglected in the PY approximation? We regard 
this as unlikely, since the two approximations are quite similar in their handling of the 
correlation functions of the fluid, giving the same functional form for c(r) in the sticky 



limit 



42| . What we deem most probable, instead, is that the failure is a combined effect 



of the approximate closure (j7|) and the consistency requirement (jSJ). Indeed, the closure is 
likely to miss a number of singular features, i.e., discontinuities and Dirac-5 peaks, that show 
up in the real correlation function js), |4l| : the isothermal compressibility, being dependent 
on the volume integral of the same function, suffers from the same defect. On the contrary, 
the internal energy of the fluid depends on the contact value of the cavity function alone: 
when thermodynamic consistency is imposed, the latter has to compensate for the missing 
contributions, unnaturally pushing away the liquid-vapor phase boundary. We stress that 
at the moment this is a mere, untested hypothesis. 

One may hope that the strong dependence of the results on the high-density boundary 
condition found for narrow SW potentials will become less drastic turning to a different form 
of the interaction. In fact, as observed above, for the HCY potential this unwelcome feature, 
albeit still there, is nevertheless much less pronounced. If this is actually the case, then the 
fully numerical solution of the SCOZA which has been developed here and the one previously 



developed in 



18| , could still be usefully applied to different forms of short-r ang e potentials, 



not to mention the instance of interactions that are not of short range 47j. Moreover, 



this solution could be he 
reference theory (HRT) 



pful in the context of the renormalization-group based hierarchical 



481 ] . specifically for enforcing the core condition on g(r) exactly (of 



course within the limit of numerical accuracy). In fact, in the smooth cut-off formulation of 



the HRT 



49j | . similarly to the SCOZA, this condition can be (and has been) implemented 



analytically for a HCY potential, but for a generic tail potential a numerical solution would 
again be needed. Also, this is always the case for the sharp-cut off HRT, irrespective of 
the form of the tail potential. In the latter theory, the core condition was implemented 
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by introducing a number of approximations 50(, which however affect the accuracy of this 



131 ] . A better procedure such as that developed here 



approach for short-range interactions 
would certainly be more adequate. 

As for the sticky limit, the results presented in this work suggest that, if thermodynamic 
consistency is to be employed as a way to improve upon the accuracy of the PY solution for 
the SHS potential ([2]), one should first get rid of the approximation common to the SCOZA 
and the PY equation: namely, that the contribution to the direct correlation function c(r) 
due to the tail interaction w(r) has the same range as w(r) itself, and therefore vanishes 
wherever w(r) does. In this respect, it is worthwhile recalling that the same approximation 
is used in conjunction with thermodynamic consistency not only in the SCOZA, but also 



in the HRT 



49| | . Therefore, if the combination of these ingredients is bound to become 



dangerous for vanishing attraction range, it is likely that this would affect the HRT as well, 
and that both the SCOZA and the HRT would benefit from a solution to this problem. 
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Appendix A: Finite-difference discretization of the SCOZA PDE 



Equation ( TTTT) is a nonlinear, diffusive-like PDE, where the quantity D(p,u) given by 
Eq. (1T21) plays the role of the reciprocal of the diffusion coefficient. This PDE is discretized 
in the following way: 



n+l/2x U j U j _ Pj 



A/3 



(A P y 



and 



Xred 



n+l 



+ 

-)' 



(A P y 



(Al) 



(A2) 
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where subscripts refer to the density p and superscripts to the inverse temperature (5 in 
the (p, 0) grid, so for instance = u(pj,/3 n ), and A(3, Ap are respectively the inverse 
temperature and density steps used in the discretization. 

Equation flAlj) has been solved by an implicit algorithm as follows: assuming that Uj is 
known at a certain inverse temperature (3 n for every pj, a guess for at the "new" tem- 
perature Pn+i is given. This is determined by solving Eqs. (jSJ), © with the approximation 
K n+ i ~ K n + A/3. The values of the inverse reduced compressibility (1/Xrcd)j +1 correspond- 
ing to are then obtained by solving Eqs. ([6]), (J7J) at fixed density p and internal energy 
u via the method described in Sec. [Til and Appendix [Bl These values are used in Eq. (|A2[) to 
obtain a trial value of D(pj,u™ +1 ^ 2 ), which is inserted into Eq. (lAip to give a closed linear 
system of equations for at /3 n+ i. The values of thus obtained replace the initial 
guess, and the procedure is iterated until convergence is obtained. The u" +1 at convergence 
are fed back into the algorithm to obtain w™ +2 , and so on. 



Appendix B: Numerical optimization of the free energy functional 

The most straightforward method to minimize the functional S defined by Eq. (1231) is 
the steepest descent 26|], in which <p{r) and K are updated iteratively by moving downhill 
in the direction opposite to that of the gradient of S, namely 

<f>i+i{r) = <pi(r) -\Ahi(r) r < a , (Bl) 

k; +1 = k; -XAv t . (B2) 

Here the indexes I, I + 1 refer to the values of <fr(r) and K at two consecutive iterations, A 
(not to be confused with the state-dependent function A(r/, r) introduced in Sec. HVl) is a 
parameter determining the step size, and we have set 

AV = -pJd 3 rAh(r)w*(r) + AU* . (B3) 

The asterisks in Eqs. (jB2|) . flB3|) denote reduced quantities K* = Ke, AU* = AU/e, w*(r) = 
w(r)/e, where e is the energy scale of w(r). For the SW potential considered here, e coincides 
with the well depth. Note that, for A to be the same dimensionless quantity in Eqs. (IB II) 
and ( 1B2I) . the differentiation of S has to be performed with respect to K* and v /p0(r). 

To implement the algorithm, one starts from a guess for K and 0(r) for r < er, and 
determines <f>(r) for r > a by the relation 0(r) = —Kw(r). In the numerical solution 
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of the SCOZA, the guess at the inverse temperature fi n +i and density pj was given by 
<f>(r,pj,/3 n +i) - <f>(r,Pj,/3n), K(pj,f3 n+1 ) ~ K(p j ,/3 n ) + A/3. By switching to reciprocal 
space, the Fourier transform of Ah(r) is determined via the OZ equation as 



Ah{k) = S HS (k)S(k)^(k) 



(B4) 



where S(k) is given by Eq. ( TT9"j) . Equation (1B4P is then inverse-transformed to give Ah(r). 
This in turn is used to determine AV by Eqs. (jTOl) . (I2TI) . (1B3|) . and the A/i(r) and AV^ thus 
obtained are used in Eqs. fIBll) . flB2j) to update 0(r) and K iteratively, until convergence is 
achieved. 

As observed in 26j, a valuable improvement over the steepest descent is the conjugate 



gradient algorithm, which we have adopted here. In this method, the direction of descent 
at the step / + 1 does not coincide with that of the gradient of S at step /, but is instead 
given by a linear combination of the gradients at all previous steps according to 



<f>i+i(r) = <M r ) - x li(r) 



k; +1 = Kt - a r, 



r < a 



(B5) 
(B6) 



where 7z(r) and Ti are defined by recurrence as 26] 



7;(r) = Ahiir) + aai^iir) 
T l = AVi + air,-! , 



r < a 



(B7) 
(B8) 



with «o = and a, for / > 1 given by 



Oil 



p [ dhAhtir) [Ahi(r) - Ah^(r)] + AV t (A^ - AV^) 

J r<l 



(B9) 



p d 3 r [Ah^r)]* + (W-i] 



r<l 



The numerator and denominator of Eq. (1B9I) are scalar products of the form {vi\vi — i>j_i), 
(vi-i\vi-i) for the present case in which the "vector" v\ consists of both the function 
y/pAhi(r) and the scalar AV. 

In the implementation of the algorithm, special attention has been paid to the choice of 
the parameter A. In both the steepest descent and conjuga te g radient methods, the optimal 
value of A at each step is determined by line minimization 26|. In the case of the conjugate 



gradient, this amounts to minimizing S[4>i(r) — X^i(r),Kl — XTi] as a function of A for 
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fixed 0;(r), ji(r), Kf, Ti. A possible way would be to use some method based on repeated 
evaluation of the functional S along the line 4>i(r) — Xj^r), — AIY We found that this 
procedure is not recommended here for a number of reasons: first, it is not very efficient, since 
it requires many evaluation of S at each step. Second, at low temperature (see Sec. [TT]) there 
is a density interval where the argument of the logarithm in Eq. ( Ilfip is not positive definite 
as a function of k, so great care has to be taken to prevent the functional S from being 



evaluated for these states 



26j . Finally, there is a subtle, yet important reason related to the 



fact that, for the hard-core plus tail potential of Eq. ([5]), both c(r) and h(r) are discontinuous 
at r = a, and display additional discontinuities if w (r) is itself discontinuous, as in the case 
of the SW potential. In order to determine with good accuracy the Fourier transform of 
the correlation functions needed in the iterative procedure, it is then advisable to subtract 
off the discontinuous contribution before performing the numerical Fourier transform, and 
add to the result its Fourier transform determined analytically, as illustrated in Sec. HT1 
However, when the functional /rung of Eq. (fT6|) is discretized in the numerical calculation, the 
fundamental relations (1T71) . (1T81) will hold rigorously if the numerical Fourier transforms of 
0(r) and Ah(r) are performed on the full functions, rather than just on their continuous part. 
As a consequence, the gain in accuracy achieved by splitting the functions to be transformed 
entails as a drawback that the solution of the Euler-Lagrange equation Sfm ng /S4>(r) = for 
r < a will not rigorously be anymore a minimum of /rung as defined by Eq. f|T6|) . and it 
is not obvious to us whether Eq. ( [TBI can be modified so as to recover this property. This 
is one more reason why we chose to perform the minimization with respect to A by an 
algorithm which would not rely on the direct evaluation of S, and hence of /king- Instead, 
we set ^(A) = S[(pi(r) — X^fi(r), K* — AI^], and determined the line minimum by solving the 
equation dip/dX = with respect to A. This was done by the Raphson- Newton method, 
namely starting from a guess value Aq and updating A recursively by the sequence 



dtp 




dX 


A fe 


(Pip 




dX 2 


A,, 



(BIO) 



In the present case we have 



d x / r w 7(r) + 57F r ' (B11) 
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§ - M"wj$m*w + vfrm-'m-M + ^ v2 ■ (B12 » 

where we have understood that 7(7") is set identically to zero for r > a, and for brevity we 
have omitted the index I in 4>(r), j(r), K*, T. All the derivatives of S are evaluated at 
M r ) - hli(r), K* - \ k T h and are easily calculated by Eqs. (ITTjl , (fTS]) . (120]) , (122]) , (123} . The 
quantities required are those already needed in Eq. (1B4I) for Ah(k). Specifically, Eqs. ( IBllj) . 
flBT2|) become 

-^t=p[d 3 rAh(rh(r) + TAV, (B13) 

^ = P Swf ^ k) ^ 2 " 2pT S^f S 2 (^*W7W + pry'^[5(A;K(A ; )] 2 (B14) 

In practice, we have found that, for each conjugate gradient step flB5j) . (1B6I) . (1B7I) . (1B8|) . a 
single Raphson-Newton step (IBlOp . (1B13I) . (1B14|) for A is sufficient. Iterating Eq. (IBlOj) so 
as to get closer to the minimum of did not lead to any significant improvement in the 
overall speed of convergence of the algorithm. 

We notice that at the critical point and on the spinodal curve, where the compressibility 
Xred diverges, the integrals that appear in Eq. flB14j) diverge as well, because S 2 (k) develops 
a non integrable singularity at k — 0, see Eq. (fT9l) . In order to account for this behavior, 
for large values of Xred the contribution to the integrals in the small-/c domain is determined 
analytically by replacing S(k) with its Ornstein-Zernike form S(k) ~ l/(x^d + ^ 2 )' anc ^ 
evaluating the remainder of the integrand at k = 0. 



Appendix C: Numerical calculation of the hard-sphere correlation function at high 
density 

We start from the exact relation 

h(r) = c(r) + J ^ ^ [h(k) - c(k)} = c{r) + ^ jf dk ^(At)^-^ , (CI) 

where c(k) has been subtracted off from h(k) in order to make the integrand decrease more 
rapidly at large k. The function ((k) = kc^ s (k)/[l — pdfts(k)], which is known analytically 
for the Waisman parametrization, is then sampled with a step Ak ~ 5 x 10~ 5 a _1 , much 
smaller than that used in the numerical Fourier transform, so that the "pseudo-Bragg" 
peaks can be included in the sampling. Each peak such that the structure factor exceeds 
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the threshold value Sns(k) = 3 is modeled by a Lorentzian A/[l + a 2 (k — ko) 2 }, where the 
parameters ko, A, a are determined so as to reproduce the position of the peak, its height, 
and its curvature. The sum of the Lorentzians is subtracted off C{k), resulting in a smooth 
function that can be transformed numerically without substantial loss of information, while 
the Fourier transform of the Lorentzians is determined analytically. In the present case in 
which ko ^ this cannot be done exactly, but we found that for narrow peaks such that 
a ^> 1 like those considered here, a satisfactory approximation is given by 

dk- S1 9 n /, fcr) , ~ J— S in(k r)e- r / a (C2) 



27r 2 r J 1 + a 2 (k — k ) 2 2-irar 



1 



7T 

cos(k r)Ci(kor) — sin(/c r) — — Si(/co r ) 



2n 2 a 2 r~~ v ""' ' " v " v " ' ~~" v ~ u " ' L2 
where Si(x) and Ci(x) are respectively the sine and cosine integrals defined as: 

Si(x) = jfdt^S, (C3) 

Ci(x) = fdt CQS(t) - 1 + ln(x) + 7 , (C4) 
Jo t 

where 7 = 0.5772 ... is Euler constant. 
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